# include "scene.h"

#ifdef SOLVE_ENERGY
# include <stdio.h>
# include <stdlib.h>

# ifdef _OMP_TEMP_BC_
# include <omp.h>
# endif

extern short interpol_1dim(const double *x, const double *u, unsigned short n, double xp, double *up);
extern double interpol_2dl(const double *x, const double *y, const double *u, double xp, double yp);
extern double thm_cond(double T);

extern void temp_bc(const struct grid *g, double ***T, short tflag, unsigned short iRK)
   {
    double *r, *z, rr[4], zz[4], (*TT)[3], ***T0, (*_T)[TEMP_INTRP_N], _z[TEMP_INTRP_N];
    short i, j, k, nr, nth, nz, _i, _j, tid, nTRD;

#   ifdef TRACE_FILE
    fprintf(g->trptr, "Entering into 'temp_bc()'...\n");
#   endif

   if( (iRK > 0 && tflag != 0) || (iRK == 0 && tflag != 1) || (iRK > 3) )
       {
        fprintf(g->lptr, "\nERROR: In 'temp_bc()': (tflag, iRK) = (%d, %d).\nThis function is designed to handle only following cases:\n\t[tflag = 0, iRK = 1, 2, 3] and [tflag = 1, iRK = 0]\n", tflag, iRK);
        exit(-1);
       }

    nr = g->nr;  nth = g->nth;  nz = g->nz; r = g->r; z = g->z; nTRD = 1;

    nTRD = 1;
#   ifdef _OMP_TEMP_BC_
    nTRD = g->nTRD;
#   endif

    _T = (double (*)[TEMP_INTRP_N])malloc(nTRD*TEMP_INTRP_N*sizeof(double));
    TT = (double (*)[3])malloc(nTRD*3*sizeof(double));

    for(_z[0] = z[NzG], k = 1; k < TEMP_INTRP_N; k++) _z[k] = 0.5*(z[NzG+k-1]+z[NzG+k]);

    T0 = g->T;
    tid = 0;
#   ifdef _OMP_TEMP_BC_
    #pragma omp parallel private(_i,_j,i,j,k,tid)
#   endif
     {
   /* Dirichlet BC at the INLET. */
#     ifdef _OMP_TEMP_BC_
      tid = omp_get_thread_num();
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < NrG+g->nr_nozzle; i++) for(j = NthG; j < nth+NthG; j++)
         {
          _T[tid][0] = T_INLET;
          for(k = 1; k < TEMP_INTRP_N; k++) _T[tid][k] = T[i][j][NzG+k-1];

          interpol_1dim(_z, _T[tid], TEMP_INTRP_N, 0.5*(z[NzG]+z[NzG-1]), &T[i][j][NzG-1]);
         }

   /* At the SOLID BASE: Constant heat-flux, SOLID_BASE_HF: Numan BC. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG+g->nr_nozzle; i < NrG+nr-g->nr_gap; i++) for(j = NthG; j < nth+NthG; j++)
         T[i][j][NzG-1] = T[i][j][NzG] - SOLID_BASE_HF*0.5*(z[NzG+1]-z[NzG-1])/thm_cond(0.5*(T0[i][j][NzG-1]+T0[i][j][NzG]));

   /* At OUTLET: Orlansky BC. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(_i = 0; _i < g->nr_gap; _i++) for(i = _i+NrG+nr-g->nr_gap, _j = 0, j = NthG; _j < nth; j++, _j++)
      switch(iRK)
           {
            case 0:
            T[i][j][NzG-1] = g->outBC.tnext_Tg[_i][_j];
            break;

            case 1:
            T[i][j][NzG-1] = g->outBC.T[2][_i][_j][0]/12 - g->outBC.T[1][_i][_j][0]/2 + 5*g->outBC.T[0][_i][_j][0]/4\
                            + g->outBC.tnext_Tg[_i][_j]/6;
            break;

            case 2:
            T[i][j][NzG-1] = g->outBC.T[2][_i][_j][0]/6 - 3*g->outBC.T[1][_i][_j][0]/4 + g->outBC.T[0][_i][_j][0]\
                            + 7*g->outBC.tnext_Tg[_i][_j]/12;
            break;

            case 3:
            T[i][j][NzG-1] = g->outBC.T[2][_i][_j][0]/24 - g->outBC.T[1][_i][_j][0]/8 + 5*g->outBC.T[0][_i][_j][0]/8\
                            + 11*g->outBC.tnext_Tg[_i][_j]/24;
            break;
           }

   /* At the TOP-SURFACE: A constant heat-flux 'TOP_HEATER_HF' in the central heater-zone and TOP_ANNULUS_HF in
      rest of the annulus-zone. */
      k = NzG+nz;
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < nth+NthG; j++)
         T[i][j][k] = T[i][j][k-1] + (i < NrG+g->nr_heater ? (g->t >= t_HEATER_ON ? TOP_HEATER_HF:0):TOP_ANNULUS_HF)\
                            *0.5*(z[k+1]-z[k-1])/thm_cond(0.5*(T0[i][j][k-1]+T0[i][j][k]));

   /* At the SIDE-WALL: A constant heat-flux "SIDE_WALL_HF". */
      i = NrG+nr;
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided)
#     endif
      for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < nth+NthG; j++)
        T[i][j][k] = T[i-1][j][k] - SIDE_WALL_HF*0.5*(r[i+1]-r[i-1])/thm_cond(0.5*(T0[i-1][j][k]+T0[i][j][k]));

   /* For ghost cells having negative radius: applying symmetric B.C. across the center.*/
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided)
#     endif
      for(k = 0; k < nz+2*NzG; k++) for(j = NthG; j < NthG+nth; j++) for(i = 0; i < NrG; i++)
         T[NrG-i-1][j][k] = T[NrG+i][j][k];

   /* For ghost cells in theta-direction: Applying symmetric B.C. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
         T[i][NthG-1-j][k] = T[i][NthG+j][k];
         T[i][j+NthG+nth][k] = T[i][NthG+nth-1-j][k];
        }

   /* Interpolating at the lower-corner cells. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp single
#     endif
        {
         rr[0] = rr[2] = 0.5*(r[NrG+nr-1]+r[NrG+nr]);
         rr[1] = rr[3] = 0.5*(r[NrG+nr]+r[NrG+nr+1]);
         zz[0] = zz[1] = 0.5*(z[NzG]+z[NzG+1]);
         zz[2] = zz[3] = 0.5*(z[NzG-1]+z[NzG]);
        }

#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided)
#     endif
      for(j = NthG; j < NthG+nth; j++)
        {
         TT[tid][0] = T[NrG+nr-1][j][NzG];
         TT[tid][1] = T[NrG+nr][j][NzG];
         TT[tid][2] = T[NrG+nr-1][j][NzG-1];
         T[NrG+nr][j][NzG-1] = interpol_2dl(rr, zz, TT[tid], rr[3], zz[3]);
        }

   /* Interpolating at the upper-corner cells. */
#     ifdef _OMP_TEMP_BC_
      #pragma omp single
#     endif
        {
         zz[0] = zz[3] = 0.5*(z[NzG+nz]+z[NzG+nz+1]);
         zz[1] = zz[2] = 0.5*(z[NzG+nz-1]+z[NzG+nz]);
        }

#     ifdef _OMP_TEMP_BC_
      #pragma omp for schedule(guided)
#     endif
      for(j = NthG; j < NthG+nth; j++)
        {
         TT[tid][0] = T[NrG+nr-1][j][NzG+nz];
         TT[tid][1] = T[NrG+nr-1][j][NzG+nz-1];
         TT[tid][2] = T[NrG+nr][j][NzG+nz-1];
         T[NrG+nr][j][NzG+nz] = interpol_2dl(rr, zz, TT[tid], rr[3], zz[3]);
        }
     }

    free(_T); free(TT);

#   ifdef TRACE_FILE
    fprintf(g->trptr, "...'temp_bc()' ends.\n");
#   endif

    return;
   }

#endif
